# Replication File for
# "Behavioral Consequences of Open Candidate Recruitment"
# Authors: Jochen Rehmert

### set directory, read data, load packages
setwd("...")

# packages and functions needed
library(foreign);library(ordinal);library(MASS);library(stargazer)
char = function(x){as.character(x)} 
num = function(x){as.numeric(char(x))}

# read data
dat <- read.dta("returnData.dta")

# subset data
dat <- dat[dat$come_home != "",]
dat$come_home <- factor(dat$come_home, levels = c("never","once a month","two three times a month",
                                                    "once a week","every other day","every day"))
# prepare independent variables
dat$LDP <- as.numeric(dat$party_id == "LDP")
dat$DPJ <- as.numeric(dat$party_id == "DPJ")
dat$SDP <- as.numeric(dat$party_id == "SDP")
dat$JRP <- as.numeric(dat$party_id == "JRP")
dat$YP <- as.numeric(dat$party_id == "YP")
dat$NPD <- as.numeric(dat$party_id == "NPD")
dat$previous <- as.numeric(dat$incumbency == "previous")
dat$freshman <- as.numeric(dat$incumbency == "freshman")
dat$ln_distance <- log(dat$distance_diet+ 1)

# kobo selected first-timer
dat$kobo_f <- NA
dat$kobo_f[dat$kobo == 1 & dat$term == 0] <- 1
dat$kobo_f[dat$kobo == 1 & dat$term >= 1] <- 0
dat$kobo_f[dat$kobo ==  0] <- 0
# kobo selected non-first-timer
dat$kobo_nf <- NA
dat$kobo_nf[dat$kobo == 1 & dat$term >= 1] <- 1
dat$kobo_nf[dat$kobo == 1 & dat$term == 0] <- 0
dat$kobo_nf[dat$kobo ==  0] <- 0

# interaction terms
dat$distKobo = dat$kobo * dat$distbirth
dat$distKobo_f = dat$kobo_f * dat$distbirth
dat$prefKobo = dat$kobo * dat$prefbirth
dat$prefKobo_f = dat$kobo_f * dat$prefbirth
dat$termKobo = dat$term * dat$kobo

kobo.party <- c("DPJ","LDP","JRP","YP", "NPD")


###############################################################################
#                              Table 1                                        #
###############################################################################

dat$home_threefold <- ""
dat$home_threefold[dat$come_home == "never" | dat$come_home == "once a month" | dat$come_home == "two three times a month"] <- "less than once a week"
dat$home_threefold[dat$come_home == "once a week"] <- "once a week"
dat$home_threefold[dat$come_home == "every other day" | dat$come_home == "every day"] <- "more than once a week"

table(dat$home_threefold, dat$come_home)

library(nnet)
#set baseline category
dat$home_threefold <- factor(dat$home_threefold, levels = c("less than once a week","once a week", "more than once a week"))
dat$home_threefold <- relevel(dat$home_threefold, ref = "once a week")

summary(third.1 <- multinom(home_threefold ~ ln_distance + party_office +  executive_office + prefbirth + margin + pty_voteshare_t0 +   kobo_f +  + LDP + YP, data  = dat[dat$party_id %in% kobo.party,]))
summary(third.2 <- multinom(home_threefold ~ ln_distance + party_office +  executive_office + prefbirth + margin + pty_voteshare_t0 +   term + kobo_f  + LDP + YP, data  = dat[dat$party_id %in% kobo.party,]))
summary(third.3 <- multinom(home_threefold ~ ln_distance + party_office +  executive_office + prefbirth + margin + pty_voteshare_t0 +   term + kobo +   LDP + YP, data  = dat[dat$party_id %in% kobo.party,]))
summary(third.4 <- multinom(home_threefold ~ ln_distance + party_office +  executive_office + prefbirth + margin + pty_voteshare_t0 +   term + kobo + termKobo +  LDP + YP, data  = dat[dat$party_id %in% kobo.party,]))

# output for table
stargazer(third.1, third.2, third.3, third.4, omit = c("LDP","YP"))

# conditional standard error for interaction term
sqrt(  diag(vcov(third.4))[8] + diag(vcov(third.4))[12] + vcov(third.4)[12,8]  ) # less than once a week
sqrt(  diag(vcov(third.4))[20] + diag(vcov(third.4))[24] + vcov(third.4)[24,20]  ) # more than once a week

###############################################################################
#                           Figure 1a and 1b                                  #
###############################################################################
dat.model <- dat[which(complete.cases(dat[,c("home_threefold","ln_distance", "prefbirth")])),c("home_threefold","ln_distance","party_office","executive_office","term","kobo", "kobo_f","margin","prefbirth","termKobo","prefKobo", "pty_voteshare_t0", "LDP","JRP","YP","NPD")]

reps = 1000
seed = 1104
x.prime = "kobo"
cond.var = "term"
int.var = "termKobo"
require(nnet)
set.seed(seed)
qoi.0.week <- c();qoi.0.less <- c()
qoi.0.more <- c();qoi.1.week <- c()
qoi.1.less <- c();qoi.1.more <- c()
qoi.2.week <- c();qoi.2.less <- c()
qoi.2.more <- c();ses <- c()
dat.tmp <- dat.model
form <- third.4[["call"]][[2]]
for(ii in 1:reps){
  print(ii)
  #data.rs <- dat.tmp[-ii,]
  data.rs <- dat.tmp[sample(1:nrow(dat.tmp), nrow(dat.tmp), replace = TRUE),]
  mod.fit <- multinom(form, data = data.rs)
  # get standard errors
  #tmp.se <- sqrt(diag(vcov(mod.fit, model = "zero")))
  tmp.se <- summary(mod.fit)$coefficients[,2]
  ses <- rbind(ses, tmp.se)
  
  # get model frame and prepare variable of interest
  tmp.dat.0.0 <- data.rs
  tmp.dat.0.0[[x.prime]] <- 0
  tmp.dat.0.0[[cond.var]] <- 0
  tmp.dat.0.0[[int.var]] <- 0 * 0
  
  tmp.dat.1.0 <- data.rs
  tmp.dat.1.0[[x.prime]] <- 1
  tmp.dat.1.0[[cond.var]] <- 0
  tmp.dat.1.0[[int.var]] <- 0 * 1
  
  tmp.dat.1.1 <- data.rs
  tmp.dat.1.1[[x.prime]] <- 1
  tmp.dat.1.1[[cond.var]] <- 1
  tmp.dat.1.1[[int.var]] <- 1 * 1
  
  tmp.dat.0.1 <- data.rs
  tmp.dat.0.1[[x.prime]] <- 0
  tmp.dat.0.1[[cond.var]] <- 1
  tmp.dat.0.1[[int.var]] <- 1 * 0
  
  tmp.dat.1.2 <- data.rs
  tmp.dat.1.2[[x.prime]] <- 1
  tmp.dat.1.2[[cond.var]] <- 2
  tmp.dat.1.2[[int.var]] <- 2 * 1
  
  tmp.dat.0.2 <- data.rs
  tmp.dat.0.2[[x.prime]] <- 0
  tmp.dat.0.2[[cond.var]] <- 2
  tmp.dat.0.2[[int.var]] <- 2 * 0
  
  
  preds00 <- data.frame(predict(mod.fit, newdata = tmp.dat.0.0[,-1], type ="prob"))
  preds10 <- data.frame(predict(mod.fit, newdata = tmp.dat.1.0[,-1], type ="prob"))
  preds11 <- data.frame(predict(mod.fit, newdata = tmp.dat.1.1[,-1], type ="prob"))
  preds01 <- data.frame(predict(mod.fit, newdata = tmp.dat.0.1[,-1], type ="prob"))
  preds02 <- data.frame(predict(mod.fit, newdata = tmp.dat.0.2[,-1], type ="prob"))
  preds12 <- data.frame(predict(mod.fit, newdata = tmp.dat.1.2[,-1], type ="prob"))
  
  week.1.0 <- mean(preds10[,1], na.rm = TRUE)
  less.1.0 <- mean(preds10[,2], na.rm = TRUE)
  more.1.0 <- mean(preds10[,3], na.rm = TRUE)
  week.1.1 <- mean(preds11[,1], na.rm = TRUE)
  less.1.1 <- mean(preds11[,2], na.rm = TRUE)
  more.1.1 <- mean(preds11[,3], na.rm = TRUE)
  week.1.2 <- mean(preds12[,1], na.rm = TRUE)
  less.1.2 <- mean(preds12[,2], na.rm = TRUE)
  more.1.2 <- mean(preds12[,3], na.rm = TRUE)
  week.0.0 <- mean(preds00[,1], na.rm = TRUE)
  less.0.0 <- mean(preds00[,2], na.rm = TRUE)
  more.0.0 <- mean(preds00[,3], na.rm = TRUE)
  week.0.1 <- mean(preds01[,1], na.rm = TRUE)
  less.0.1 <- mean(preds01[,2], na.rm = TRUE)
  more.0.1 <- mean(preds01[,3], na.rm = TRUE)
  week.0.2 <- mean(preds02[,1], na.rm = TRUE)
  less.0.2 <- mean(preds02[,2], na.rm = TRUE)
  more.0.2 <- mean(preds02[,3], na.rm = TRUE)
  
  ame.week.0 <- week.1.0 - week.0.0
  ame.more.0 <- more.1.0 - more.0.0
  ame.less.0 <- less.1.0 - less.0.0
  
  ame.week.1 <- week.1.1 - week.0.1
  ame.more.1 <- more.1.1 - more.0.1
  ame.less.1 <- less.1.1 - less.0.1
  
  ame.week.2 <- week.1.2 - week.0.2
  ame.more.2 <- more.1.2 - more.0.2
  ame.less.2 <- less.1.2 - less.0.2
  
  qoi.0.less <- rbind(qoi.0.less, ame.less.0)
  qoi.0.more <- rbind(qoi.0.more, ame.more.0)    
  qoi.0.week <- rbind(qoi.0.week, ame.week.0)
  qoi.1.less <- rbind(qoi.1.less, ame.less.1)
  qoi.1.more <- rbind(qoi.1.more, ame.more.1)      
  qoi.1.week <- rbind(qoi.1.week, ame.week.1)
  qoi.2.less <- rbind(qoi.2.less, ame.less.2)
  qoi.2.more <- rbind(qoi.2.more, ame.more.2)     
  qoi.2.week <- rbind(qoi.2.week, ame.week.2)
}

third.4.boots <- list(qoi_0_less = qoi.0.less, qoi_0_week = qoi.0.week, qoi_0_more = qoi.0.more,
                      qoi_1_less = qoi.1.less, qoi_1_week = qoi.1.week, qoi_1_more = qoi.1.more,
                      qoi_2_less = qoi.2.less, qoi_2_week = qoi.2.week, qoi_2_more = qoi.2.more)



ggData <- data.frame(rbind(cbind(ame = mean(third.4.boots$qoi_0_less, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_0_less, probs = 0.025), upper = quantile(third.4.boots$qoi_0_less, probs = 0.975),  level = "less than once a week", term = 0),
                           cbind(ame = mean(third.4.boots$qoi_0_week, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_0_week, probs = 0.025), upper = quantile(third.4.boots$qoi_0_week, probs = 0.975), level = "once a week", term = 0),
                           cbind(ame = mean(third.4.boots$qoi_0_more, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_0_more, probs = 0.025), upper = quantile(third.4.boots$qoi_0_more, probs = 0.975),level = "more than once a week", term = 0),
                           cbind(ame = mean(third.4.boots$qoi_1_less, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_1_less, probs = 0.025), upper = quantile(third.4.boots$qoi_1_less, probs = 0.975), level = "less than once a week", term = 1),
                           cbind(ame = mean(third.4.boots$qoi_1_week, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_1_week, probs = 0.025), upper = quantile(third.4.boots$qoi_1_week, probs = 0.975), level = "once a week", term = 1),
                           cbind(ame = mean(third.4.boots$qoi_1_more, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_1_more, probs = 0.025), upper = quantile(third.4.boots$qoi_1_more, probs = 0.975), level = "more than once a week", term = 1),
                           cbind(ame = mean(third.4.boots$qoi_2_less, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_2_less, probs = 0.025), upper = quantile(third.4.boots$qoi_2_less, probs = 0.975), level = "less than once a week", term = 2),
                           cbind(ame = mean(third.4.boots$qoi_2_week, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_2_week, probs = 0.025), upper = quantile(third.4.boots$qoi_2_week, probs = 0.975), level = "once a week", term = 2),
                           cbind(ame = mean(third.4.boots$qoi_2_more, na.rm = TRUE), ci = 95, lower = quantile(third.4.boots$qoi_2_more, probs = 0.025), upper = quantile(third.4.boots$qoi_2_more, probs = 0.975), level = "more than once a week", term = 2),
                           cbind(ame = mean(third.4.boots$qoi_0_less, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_0_less, probs = 0.05), upper = quantile(third.4.boots$qoi_0_less, probs = 0.95),  level = "less than once a week", term = 0),
                           cbind(ame = mean(third.4.boots$qoi_0_week, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_0_week, probs = 0.05), upper = quantile(third.4.boots$qoi_0_week, probs = 0.95), level = "once a week", term = 0),
                           cbind(ame = mean(third.4.boots$qoi_0_more, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_0_more, probs = 0.05), upper = quantile(third.4.boots$qoi_0_more, probs = 0.95),level = "more than once a week", term = 0),
                           cbind(ame = mean(third.4.boots$qoi_1_less, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_1_less, probs = 0.05), upper = quantile(third.4.boots$qoi_1_less, probs = 0.95), level = "less than once a week", term = 1),
                           cbind(ame = mean(third.4.boots$qoi_1_week, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_1_week, probs = 0.05), upper = quantile(third.4.boots$qoi_1_week, probs = 0.95), level = "once a week", term = 1),
                           cbind(ame = mean(third.4.boots$qoi_1_more, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_1_more, probs = 0.05), upper = quantile(third.4.boots$qoi_1_more, probs = 0.95), level = "more than once a week", term = 1),
                           cbind(ame = mean(third.4.boots$qoi_2_less, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_2_less, probs = 0.05), upper = quantile(third.4.boots$qoi_2_less, probs = 0.95), level = "less than once a week", term = 2),
                           cbind(ame = mean(third.4.boots$qoi_2_week, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_2_week, probs = 0.05), upper = quantile(third.4.boots$qoi_2_week, probs = 0.95), level = "once a week", term = 2),
                           cbind(ame = mean(third.4.boots$qoi_2_more, na.rm = TRUE), ci = 90, lower = quantile(third.4.boots$qoi_2_more, probs = 0.05), upper = quantile(third.4.boots$qoi_2_more, probs = 0.95), level = "more than once a week", term = 2)))

colnames(ggData)[1] <- "ame"
ggData[,1] <- num(ggData[,1])
ggData[,2] <- num(ggData[,2])
ggData[,3] <- num(ggData[,3])
ggData[,4] <- num(ggData[,4])
ggData[,5] <- ordered(ggData[,5], levels = c("less than once a week","once a week", "more than once a week"))

# ----- Figure 1a ----- #

yticks <- c(-0.1,0,0.1,0.2,0.3)
xticks <- c(1.2,1.8)
par(oma = c(3,1,1,1),
    mai = c(3,1,1,1),
    mar =  c(3, 5.5, 1, 1) + 0.1,
    mgp = c(3, 1, 1), 
    xpd = NA,
    bty = 'n')
plot(ggData$ame[ggData$term == 0 & ggData$level != "once a week" & ggData$ci == 95] ~ c(1.2,1.8), 
     ylab = "Average Marginal Effect", xlab = "", pch = 19, xlim = c(1,2), ylim = c(-0.1,.4), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n')
segments(1.2, ggData$lower[ggData$term == 0 & ggData$level == "less than once a week" & ggData$ci == 95], 
         1.2,ggData$upper[ggData$term == 0 & ggData$level == "less than once a week" & ggData$ci == 95] , lty = 2   )
segments(1.8, ggData$lower[ggData$term == 0 & ggData$level == "more than once a week" & ggData$ci == 95], 
         1.8,ggData$upper[ggData$term == 0 & ggData$level == "more than once a week" & ggData$ci == 95] , lty = 2   )
segments(1.2, ggData$lower[ggData$term == 0 & ggData$level == "less than once a week" & ggData$ci == 90], 
         1.2,ggData$upper[ggData$term == 0 & ggData$level == "less than once a week" & ggData$ci == 90] , lty = 1   )
segments(1.8, ggData$lower[ggData$term == 0 & ggData$level == "more than once a week" & ggData$ci == 90], 
         1.8,ggData$upper[ggData$term == 0 & ggData$level == "more than once a week" & ggData$ci == 90] , lty = 1   )
lines(x = c(1,2), y = c(0,0), col = "red", lty = "dashed")
text(par("usr")[1],y= yticks, labels = yticks,
     pos = 2, xpd = TRUE, cex = 1.25)
text(par("usr")[3], x = xticks, labels = c("Less than\n once a week", "More than\n once a week"), 
     pos = 1, xpd = TRUE, cex = 1.25)



# ----- Figure 1b ----- #

yticks <- c(-0.1,0,0.1,0.2,0.3)
xticks <- c(1.2,1.8)
par(oma = c(3,1,1,1),
    mai = c(3,1,1,1),
    mar =  c(3, 5.5, 1, 1) + 0.1,
    mgp = c(3, 1, 1), 
    xpd = NA,
    bty = 'n')
plot(ggData$ame[ggData$term == 1 & ggData$level != "once a week" & ggData$ci == 95] ~ c(1.2,1.8), 
     ylab = "Average Marginal Effect", xlab = "", pch = 19, xlim = c(1,2), ylim = c(-0.1,.4), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n')
segments(1.2, ggData$lower[ggData$term == 1 & ggData$level == "less than once a week" & ggData$ci == 95], 
         1.2,ggData$upper[ggData$term == 1 & ggData$level == "less than once a week" & ggData$ci == 95] , lty = 2   )
segments(1.8, ggData$lower[ggData$term == 1 & ggData$level == "more than once a week" & ggData$ci == 95], 
         1.8,ggData$upper[ggData$term == 1 & ggData$level == "more than once a week" & ggData$ci == 95]  , lty = 2  )
segments(1.2, ggData$lower[ggData$term == 1 & ggData$level == "less than once a week" & ggData$ci == 90], 
         1.2,ggData$upper[ggData$term == 1 & ggData$level == "less than once a week" & ggData$ci == 90] , lty = 1   )
segments(1.8, ggData$lower[ggData$term == 1 & ggData$level == "more than once a week" & ggData$ci == 90], 
         1.8,ggData$upper[ggData$term == 1 & ggData$level == "more than once a week" & ggData$ci == 90]  , lty = 1  )
lines(x = c(1,2), y = c(0,0), col = "red", lty = "dashed")
text(par("usr")[1],y= yticks, labels = yticks,
     pos = 2, xpd = TRUE, cex = 1.25)
text(par("usr")[3], x = xticks, labels = c("Less than\n once a week", "More than\n once a week"), 
     pos = 1, xpd = TRUE, cex = 1.25)


###############################################################################
#                             Table A15                                       #
###############################################################################
summary(third.1 <- multinom(home_threefold ~ ln_distance + party_office + volatility + enc + pure_smd + executive_office + prefbirth + margin + pty_voteshare_t0 +   kobo_f +  + LDP + YP, data  = dat[dat$party_id %in% kobo.party,], model = T))
summary(third.2 <- multinom(home_threefold ~ ln_distance + party_office + volatility + enc +pure_smd+ executive_office + prefbirth + margin + pty_voteshare_t0 +   term + kobo_f  + LDP + YP, data  = dat[dat$party_id %in% kobo.party,]))
summary(third.3 <- multinom(home_threefold ~ ln_distance + party_office + volatility + enc +pure_smd+ executive_office + prefbirth + margin + pty_voteshare_t0 +   term + kobo +   LDP + YP, data  = dat[dat$party_id %in% kobo.party,]))
summary(third.4 <- multinom(home_threefold ~ ln_distance + party_office + volatility + enc +pure_smd+ executive_office + prefbirth + margin + pty_voteshare_t0 +   term + kobo + termKobo +  LDP + YP, data  = dat[dat$party_id %in% kobo.party,]))

# output for table
stargazer(third.1, third.2, third.3, third.4, omit = c("LDP","YP"))


###############################################################################
#                             Table A14                                       #
###############################################################################


summary(first.1 <- clm(come_home ~ ln_distance + party_office + executive_office + margin + term + prefbirth +kobo + pty_voteshare_t0 +  LDP + JRP + YP + NPD, data  = dat[dat$party_id %in% kobo.party,], link = "cauchit"))
summary(first.2 <- clm(come_home ~ ln_distance + party_office + executive_office + margin + termKobo +term  + kobo + prefbirth + pty_voteshare_t0 +  LDP + JRP + YP + NPD, data  = dat[dat$party_id %in% kobo.party,], link = "cauchit"))
#onditional standard error for interaction term
sqrt(diag(vcov(first.2))[10]  +  diag(vcov(first.2))[11] + 2 * vcov(first.2)[11,10])  

summary(second.1 <- clm(come_home ~ ln_distance + party_office + executive_office + margin  + prefbirth + kobo_f + pty_voteshare_t0 +  LDP + YP , data  = dat[dat$party_id %in% kobo.party,], link = "cauchit"))
summary(second.2 <- clm(come_home ~ ln_distance + party_office + executive_office + margin + term + prefbirth + kobo_f + pty_voteshare_t0 +  LDP + YP , data  = dat[dat$party_id %in% kobo.party,], link = "cauchit"))

# output table
stargazer(second.1, second.2, first.1, first.2, omit = c("LDP","YP","JRP","NPD"))


